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ABSTRACT 

Context. The vertical thickness of debris discs is often used as a measure of these systems' dynamical excitation, and as clues to the 
presence of hidden massive perturbers such as planetary embryos. However, this argument might be flawed because the observed dust 
should be naturally placed on inclined orbits by the combined effect of radiation pressure and mutual collisions. 
Aims. We critically reinvestigate this issue and numerically estimate what the "natural" vertical thickness of a collisionally evolving 
disc is, in the absence of any additional perturbing body. 

Methods. We use a deterministic collisional code, following the dynamical evolution of a population of indestructible test grains 
suffering mutual inelastic impacts. Grain differential sizes as well as the effect of radiation pressure are taken into account. 
Results. We find that, under the coupled effect of radiation pressure and collisions, grains naturally acquire inclinations of a few 
degrees. The disc is stratified with respect to grain sizes, with the smallest grains having the largest vertical dispersion and the bigger 
ones clustered closer to the midplane. 

Conclusions. Debris discs should have a minimum "natural" observed aspect ratio /i„„„ ~ 0.04±0.02 at visible to mid-IR wavelengths 
where the flux is dominated by the smallest bound grains. These values are comparable to the estimated thicknesses of many vertically 
resolved debris discs, as is illustrated with the specific example of AU Mic. For all systems with h ~ h^j,„ the presence (or absence) 
of embedded perturbing bodies cannot be inferred from the vertical dispersion of the disc. 
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' 1. Introduction 

1.1. vertical structure of debris discs 

Dusty debris discs have been observed around main sequence 
stars for over two decades. It is believed that the detected dust 
. is the lower end of a collisional cascade starting at larger, 

■ possibly planetesimal-sized bodie s which are steadily erod- 
ing due to high velocity impacts (IWvatt I l2008l) . As of April 
2009, 18 such debris discs have been spatially resolved (see 
http://www.circumstellardisks.org/), mostly in scattered light 

, (visible and near IR) but also in some cases at mid-IR and up 

■ to millimetre wavelengths. These images have revealed a great 

■ variety of radial and azimutal features, such as rings, clumps 
or two-side asymmetries, and almost no system displays a fea- 
tureless structure. These features have been extensively stud- 
ied, mostly through numerical modelling, considerably improv- 
ing our understanding of the processes at play in debris discs. 
Depending on the specificities of each system, they have been 
interpreted as being the signatures of perturbing planets or stel- 
lar companions, or the results of transitory catastrophic events, 
or due to interactions between the dust and remnant gas. 

In comparison, fewer discs have been resolved in the vertical 
direction. Edge-on seen systems are in this respect the most 
favourable cases. The best resolved disc so far is /J-Pictoris, for 
which the vertical structure is fairly well constrained, with a rel- 
atively constant thickness H ~ 15 - 20 AU (defined by half the 
vertical FWHM) in the r < 120AU region and an as pect ra- 
tioQ /i - H/r ~ 0.05 - 0.1 in the outer regions (.Heap et alj 
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120001: iGoUmowski et al ."2006'). Observations in scattered-light 
have also revealed a striking warped profile an d possibly the 
prese nce of two sUghtly tilted distinct discs (e. g lMouillet et al.l 
ll997l:lHeap etani2000l:[Goiimowski et alJ2006h . The other rela- 
tively well resolved disc is around the /? Pic coevol star Au Mic. 
It resembles to some extent a scaled-down Pic, with a roughly 
flat d isc in the inner r < 30 AU region and a H/r ~ 0.05 beyond 
that jKrist et alj|2005h . but no obvious vertical asymmetries or 
warps. For the two other resolved edge-on systems, HD15115 
and HD 139664, the vertical structure is less constrained be- 
cause of poorer resolution and signal to noise ratios. HD15115 
has a measured aspect ratio of ~ 0.05 but seems to display 
pronoun ced asymmetries in its half-width at quarter maximum 
profile ( Kalas etalj|2007h . w hile HD139 664's vertical profile 
has not been constrained yet jKalas et al.l l2006). For discs seen 
head-on or at high angle, the geometrical configuration is less 
favourable, and the vertical structure is difficult to retrieve, or 
only through model dependent inversion of brightness profiles. 
The HD181327 rings for example, might have H/r ratio as large 
as about 0. 1 at the positions of maximum surface density, but th e 
actual ratios could be two times smaller jSchneider et al.ll2006l) . 
There is however one non-edge-on system for which a precise 
thickness estimate has been provided, F omalhaut, which m ight 
have an aspect ratio as low as 0.0125 (Kalas et"ani2005h . We 
discuss this specific case in more detail in Section 4.3. 

Note that there is often an ambiguity regarding the way the 
thickness H is defined. For an edge-on seen system, this thick- 
ness, and all the directly related parameters such as aspect ratio 
or opening angle, is measured on the disc's luminosity profile 
along its two radial extensions. In other words, it is not a mea- 
sure of the disc's local thickness at a given radial location, but 
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an observed thickness, integrating along the Hne of sight grains 
from different radial distances Q. For systems seen at higher an- 
gle, however, the thickness inferred through model dependent 
fitting is usually the real local thickness at a given radial distance 
from the star 



1.2. disc thickness as a measure of dynamical excitation 

Despite the relative lack of observational data, disc vertical 
structures, or at least their global aspect ratios, have often been 
used as crucial parameters in debris disc studies, in particular 
those aimed at modelling their dynamical and collisional evolu- 
tion. Indeed, disc thickness is commonly considered as being the 
only observable that can give a direct information about the sys- 
tem's dynamical excitation, and hence about the processes and 
bodies shaping it. The argument is the following: a system of 
mutually colliding bodies orbiting around a central object tends 
toward an equilibrium where impact velocities Av are on aver- 
age of the order of the bodies surface escape velocity Vgsc- In 
addition, there is an equipartition between in-plane and verti- 
cal motions such as (/) ~ l/2{e}, where / and e are the orbital 
inclinations and eccentricities respectively. This means that if 
the system is only made of small dust grains, then it should be 
almost razor thin, with (/) ~ 0. If on the contrary it has a fi- 
nite observable thickness, then it means that "something" is dy- 
namically stirring up the system to these large values. This argu- 
ment is often used to infer the presence of unseen large bodies 
which are governing the disc's dynamics, and also to quantita- 
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taking into account several paramete rs such as the disc's age, 
surface density, etc., can be found in lOuillen et alj jlOOf). We 
only present here a simplified version, based on the first order 
assumption that t he encounter veloci ty dispersion (Av) is of the 
order of Vesd^big) d ArtvmowicJ 1 9971) . Sbig can then be estimated 
through the equipartition condition {;) ~ 1 /2(e) couple d to the 
relation dWetheriU & Stewartlll993l:IWvatt & Denll2002h : 



< Av >' 



( 



< f > +1.25 <e^ > 



j VKep(r) 



c(Sbig) 



(1) 



where vxepir) is the orbital velocity at radial distance r. For a 
typical debris disc aspect ratio Hjr - 0.05 at a typical distance 
r - 50AU from the star, this gives Vesdstig) - 400m. s"', and 
thus Sbig ~ 200 - 500 km, depending on the bodies physical den- 
sity B 

This argument gains additional credibility from the fact that 
the inferred Sbig values make sense within the frame of our un- 
derstanding of debris disc evolution. Indeed, it is believed that 
the debris disc phase starts when the bulk of the planet for- 
mation process is over and large, Ceres to Moon-sized plan - 
etary embryos have formed (e.g. iKenvon & Bromlevi |2004|) . 
Observations thus seem to corroborate theory. 

However, there could be a major problem with this 
observation-based derivation of Sbig. Indeed, it overlooks one 
crucial fact, i.e. that the smallest observable dust grains, those 
which dominate the flux in scattered light, are strongly affected 
by radiation pressure from the central star. These particles' ec- 



^ all those with r' = r.sin(8), where r is the distance measured along 
the ansae and r' the real radial distance 

^ This is a first order estimate. In reality the constraint is not directly 
on stig but on t he product Sbin^bie, where £4,5 is the surface density of 
Sbig bodies (see lOuillen et alJl2007L for more details) 



centricities are thus not solely imposed by the gravitational per- 
turbations of hypothetical large embedded bodies but also by the 
eccentricity e^p they naturally gain when being produced. For 
the simplified case of a body of size s produced from a parent 
body on a circular orbit, this eccentricity reads 



eRp - 



1 -/3(s) 



(2) 



where /3(s) ~ 0.5 x Scui/s is the ratio between the radiation pres- 
sure force and gravity (Ji > 0.5 for unbound orbits with s < Scm). 
Once placed on such high-e orbits, these particles will impact 
other bodies at high velocities. These collisions will necessar- 
ily transfer a fraction of the high horizontal (in-plane) velocities 
into vertical ones. Exactly how much Av is transfeiTed depends 
on the geometry of the impact, the mass ratio between impactor 
and target and the distribution of fragments produced by the col- 
lision, but it is certain that some in-plane to out-of -plane veloc- 
ity transfer will occur. As a consequence, the disc should have 
a natural tendency to thicken, even without any external source 
of perturbations. If proven significant, this effect could greatly 
affect the interpretations of debris disc images, by questioning 
the direct link between disc thickness and dynamical stirring. 

In this study, we propose to numerically investigate this issue 
by quantifying the "natural" vertical structure of a coUisionally 
evolving debris disc, and in particular its equilibrium vertical 
thickness. 



2. Model 

Numerical studies of debris discs fall into two main cate- 
gories: those investigating the collisional and size distribution 
evolution of the system are usually statistical particle-in-the- 
box models, with no or poor spatial resolution and dynamical 
evolution (e glKrivov et al.1 120061: iThebault & A ugereau 1 120071: 
iLohne et al.ll2008l) . while those studying the dynamics and the 
formation and evolution of spatial structures are mostly N-body 
type codes, where size distributio ns and mutual col hsions are 
usually neglected (see for example iReche et"aLll2008[) . 

The present problem poses however a specific challenge, 
inasmuch as it aims to quantify the effect of mutual colli- 
sions, within a population of bodies of different sizes, on the 
dynamical structure. Unfortunately, to study both the dynam- 
ical and the collisional (and size) evolution of a dusty disc 
is far beyond the present day computing capacities, and only 
a simplified approach is possible. In this respect, one of the 
most suitable available codes is pr obably the one develope d 3 
decades ago by iBrahid (Il977h and iBrahic & HenonI (Il977h to 
stu dy the collisional evolutio n of planetary rings, later upgraded 
by IThebault & Brahi3 ( Il998l) for the study of collisions among 
planetesimals. We chose to use here a revised version of this 
code, adapted to the present problem. 



2.1. code 

It is an N-body type deterministic code, following the dynami- 
cal evolution of a population of massless test particles under the 
influence of a central body's potential as well as possible addi- 
tional gravitational perturbers. It has a build-in collision search 
algorithm, which tracks, at each time step, all possible 2-body 
mutual encounters within the population. Due to the unavoidable 
limited number of test particles, typically a few lO'*, the statistics 
of impacts is too poor with the particles' real sizes. We thus take 
the usual "inflated radius" assumption, assigning each particle a 



Thebault: Vertical structure of debris discs 



3 



Table 1. Nominal case set up 



2.2. setup 



Radial extent^ 
Radial surface density profile 
Number of test particles A' 
Size range ^ 

Normal rebound coef. 6^ 
Transverse rebound coef. 67- 
Size distribution 
Inflated radius (in AU) 
Initial dynamical excitation^ 



10 < r < 100 AU 

E oc r-°-^ 
5x lO'' 

y3(w) = 0.025 <p{s)<l3{ 

-0.3 

1 

dN(s) oc s-^-^ds 



0.4 



" initial location of particles' periastron 

* as parameterized by y6 oc 1/^ 
the value for {edy„) is that of hypothetical large particles with /? = 
0, not taking into account the eccentricity component due to radiation 
pressure 



numerical size s„um large enough to get enough impact statistics. 
Such an approximation is justified as long as 5„„,„ is smaller than 
the minimal dista nce over which dyna mical conditions change 
in the system (e.g. lCharnoz et al1l200lh . 

Collisions are then treated as inelastic rebounds. Collision 
outcomes are estimated by computing the momentum before 
and after impact phej and pa/t in the two-bodies centre of mass 
frame. The link between p^jj, and ph^j- is parameterized by a 
normal and transversal rebound coefficients, €n and ej, such as 
En - -I and = 1 for a perfectly elastic collision. 

For the present problem, we modified the code in order to 
accommodate for a size distribution within the test particle pop- 
ulation, and more specifically for impacts between objects of dif- 
ferent sizes. In this case, the momentum balance is weighted by 
the mass ratio between the 2 impacting bodies. Another crucial 
improvement is that we now take into account the effect of radia- 
tion pressure on the smallest particles. This force is computed by 
correcting the gravitational force from the central star by a fac- 
tor (1 - /3(s)), where /5(s) = pRp/Fgrav In this respect, it might 
be more convenient to scale particle sizes by their /3 value, with 
/3{s) = 0.5 for the smallest bound objects (if released from cir- 
cular orbits). 

Of course, this "bouncing balls" model is a crude approxima- 
tion of the real behaviour of a collisionally evolving debris disc. 
In a real system high-Av collisions should result in fragmenta- 
tion rather than inelastic rebounds, producing each time a cloud 
of small fragments. As mentioned before, such chain reactions 
of fragment-producing collisions are impossible to numerically 
model. However, if we make the reasonable assumption that a 
steady-state has been reached in the system, so that the size dis- 
tribution no longer globally evolves due to collisions Q, then the 
present code gives a first order estimate of how mutual collisions 
between bodies of different sizes redistribute in and out-of-plane 
velocities. Indeed, this effect is for each impact proportional to 
the relative velocity, the mass ratio between impactors, the ge- 
ometry of the impact and energy dissipated at impact, all pa- 
rameters which our code handles (the energy dissipations being 
tuned-in by the values for eN and €t). We discuss these issues in 
more details in Sect|4.1| 



For the sake of clarity, we consider a nominal reference case, 
around which several free parameters will be explored individ- 
ually. This reference case roughly coiTesponds to a hypothetic 
typical debris disc, spatially extended from 10 to 100 AU. We 
consider N = 5 x 10"* test particles. The absolute values of their 
physical radius are not relevant here. Indeed, for the collision 
outcome procedure, the decisive factor is the impactors mass (or 
s^) ratio. For the collision search routine, the assumed inflated 
radius is proportional to the physical size s, but the proportion- 
ality factor, and thus the "real" value of s, is of no importance, 
the only constraint being on the absolute value of s„um which 
has to be large enough impact statistics and small enough not 
to introduce any bias (see below). As for the orbital evolution 
computation, what matters is the value of /3(s) for estimating 
the radiation pressure force. We thus chose to scale all particle 
sizes by their /3{s) oc 1/s value. We assume that the size distri- 
bution has reached a steady state profile following a power law 
dN oc s''ds, or more exactly dN oc /3'i^^d/3. We assume here the 
classical value q = -3.5 for idea lized infinite self similar col- 
lisional cascades (iDohnanvill 1 969 ), but another possible values 
are explored. With such steep size distributions, the size range 
which can be modelled with 5x10"* particles is necessarily lim- 
ited. We take here i,„,„ and s„,ax such as /3(s„i„) - 0.4 (close to 
the blow-out cut-off size) and /5is„ax) - 0.025. This size range 
covers however most of the crucial grain population which dom- 
inates the flux in scattered light. For the collision search routine, 
the inflated radius 5„,„„ is taken proportional to the particles real 
sizes and is set such as to give enough collision statistics without 
introducing a bias in impact rate and velocity estimates. We take 
s„u,n - lO"** AUx(J3(sii,i„) / p{s), which is a typical value for simi- 
lar im pact search routines (e.g. Charnoz et al. 2001; Xie & Zhoul 
l2008h . 

Fo r collision ou t comes , w e follow iThe bault & Brahid 
(119981) . ICharnoz et al.1 (I2OOII) and lThebault&B eust {200^ and 
take £n = -0.3 and ex - I. This corresponds to a 90% energy 
dissipation for head-on impacts, which is the usual first-order 
estimate for high v elocity impacts (e.g. lPetit & Farinellalll993t 
iFuiiwara et al.lll989 ). but other values are explored. 

In order to avoid confusion, we label (edyn) the "dynami- 
cal" component of the particles initial eccentricities, (e^/y,,) is 
thus a size independent quantifty, the other component being 
the (size dependent) radiation pressure imposed eccentricity e^p 
(see Eq|2]l. For the starting dynamical conditions, we consider 
an extreme fiducial case with (ed,-,,) = (/) = 0, i.e. the most dy- 
namically quiet state possible for the system. This is in order to 
get a minimum value for the "natural" thickening effect studied 
here, i.e., the collisonal transfer of high Av induced by radiation 
pressure. Note that even if e,/v„ = 0, small grains will have a 
non-zero initial eccentricity because of radiation pressure. 

All initial conditions are summarized in TablT] 

We let the simulations run until we do not see any further 
evolution of the average dynamical characteristics for the dif- 
ferent particle size groups. In practice this is achieved once the 
average number of impacts per particle is ~ 3. Beyond this point, 
the system has reached a steady state and only the slow and pro- 
gressive decrease of (e) and (/), due to the energy dissipation at 
impacts, is observed (see Fig|2]and Sec 14. II for a mode detailed 
discussion). In this respect, the absolute time scale is not rele- 
vant as long as the relaxation time of a few impacts per bodies 



for each given size there is an equilibrium between the particles lost 
by fragmenting collisions and the ones newly produced as fragments 
from collisions involving bigger bodies 
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Fig. 1. Final positions in the vertical plane, once the system has 
reached a collisional steady state, for the 5 x lO'* particles in 
the nominal run. The colour scale is indicative of the particles 
physical sizes. 



is short compared to the age of the system, a condition which is 
most likely to be fulfilled for any observed debris disc (see for 
instance .Wvatt ..2005i) 

3. Results 

Fig[T] shows the final positions, once the system is coUisionally 
relaxed, of all particles in the (x,z) plane. An obvious result is 
that, in sharp contrast to its initially perfectly thin structure, the 
disc is now significantly extended in the vertical direction, with a 
large fraction of the test particles being more than 10 AU above 
or below the midplane. Moreover, the disc's vertical structure 
is strongly stratified with respect to sizes. The bigger objects 
are clustered close to the midplane while the (most numerous) 
smaller grains with high fi have a much larger spread in z- This 
feature appears more clearly in Fig|2] showing the evolution of 
{i) for different size ranges: particles close to the blow-out size 
reach inclinations of about 8°, but the largest simulated grains 
(about 15 times this minimum size) settle to less than 0.5°. These 
results are logical since the smallest grains are those that have 
the most eccentric orbits because of radiation pressure. As such 
they have the highest impact velocities and thus more kinetic 
energy to transfer from the horizontal to the vertical direction. 
Likewise, since they are the smallest grains, they are the ones 
which are most likely to get large dynamical "kicks" from im- 
pacts with larger impactors (because of the momentum balance 
in the centre of frame). 

To better visualize the effect of this inclination increase on 
the disc's aspect, we consider two cases: a disc seen edge-on, 
for which we show a synthetic profile of the disc's vertical 
FWHM and FWO.IM in scattered light assuming grey scatter- 
ing (FigO, and a disc seen at higher angle, for which we display 
the azimuthally averaged local thickness (FWHM and FWO. IM) 
Hgeomif) of the disc, derived from the vertical profile of the in- 
tegrated geometrical cross section of all particles present in a 
given radial annulus (Fig|4|. Due to the limited number of par- 
ticles, the resolution is limited to 2AUxlAU and the profiles are 
relatively noisy, especially in the outermost regions because of 
the radial and vertical dilution of the particle number density. It 



Fig. 2. Evolution of {/), for different grain sizes in the nominal 
run. 



40 




-40 1 I 

-100 -50 50 100 

X (AU) 

Fig. 3. Synthetic image of a disc seen edge-on, once the steady 
collisonal state is reached, in scattered light assuming grey scat- 
tering. The dotted line corresponds to the vertical full width at 
half maximum (FWHM), while the full line marks the location 
of the full width at one tenth maximum FWO. IM 
. The fluxes have been averaged over 2AUxlAU bins. 



is however enough to resolve both the FWHM and the FWO. IM 
in the whole 0-120 AU region. As can be seen, the two different 
thicknesses, i.e., the edge-on profile derived Hedge and the local 
geometrical Hgeom have relatively comparable values, and can 
as a first approximation be considered to be identical within the 
error bars imposed by the noisiness of our profiles. The only dif- 
ferences occur close to the inner edge of the parent bodies disc 
around 10 AU, where Hge„m drops to while the edge-on seen 

luminosity has a roughly constant thickness in the whole 20 

to 20 AU region because it integrates along the line of sight con- 
tributions from particles at different radial locations. We ignore 
for the time being these limited differences (we rediscuss them 
in Sect j4.3l l and consider hereafter that the parameter H stands 
for both possible definitions of the disc's thickness. 

We see that, in both cases, the system's aspect ratio h - Hjr 
is relatively constant and is roughly ~0.04. The FWO. IM profile 
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Fig. 4. Same run as in Fig 13 but displaying the FWHM and 
FWO.IM for the azimuthally averaged local geometrical thick- 
ness as a function of radial distance 



We first consider two extreme values for normal rebound co- 
efficient €m, -0.2 and -0.5, corresponding to a 96% and 75% en- 
ergy dissipation respectively, both values being as the upper and 
lower range of expect ed energy losse s for high velocity frag- 
menting impacts ( Fuiiwara et al.|[T989 ). As could be logically 
expected, (iV) decreases with increasing energy dissipation, but 
the dependence on remains relatively weak, with only 
varying by ~ 50% between our two extreme cases. 

Surprisingly, varying the dynamical excitation of the sys- 
tem doesn't lead to any change in the final {iw) as long as 
{cdyn) ^ 0.15. For this range of (ejy,,), the inclinations acquired 
by collisional redistribution of the eccentric velocities due to ra- 
diation pressure (which do not depend on ejv«) dominate over 
the value of {idyn)- 

The only case for which we get a significantly flatter disc 
is when considering a much shallower slope for the size distri- 
bution (q - -2.5). In this case, {idyn) - 2.2°, almost half the 
inclination of the nominal case. This is because in this case the 
system's geometrical cross section, and thus the flux in scattered 
light, is now dominated by the larger bodies of the distribution 
which stay on low inclination orbits close to the mid-plane. 
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Fig. 5. Evolution of {iw), averaged over all particle sizes, for dif- 
ferent initial set-ups 

is much wider, with /iq.i ~ 0.1 -0.15. The seemingly "flatter" as- 
pect of the disc, when compared to the (x,z) positions displayed 
in FiglT] is due to the fact that in deriving the scattered light 
flux, particles' geometrical cross sections are taken into account, 
which increases the contribution of larger grains, having smaller 
z dispersion. This is illustrated in Fig|5] showing {iw), the aver- 
age inclination for the whole grain population weighted by each 
particle's cross section. It can be seen that {iw) ~ 3.5°, roughly 
half that of the smallest grains (Fig|2]i. Note that this weighted 
{iw) is approximately equal to V(2)/i, which is consistent with 
the relation between inclination and aspe ct ratio, for s mall (/) 
values, in single-sized particle discs (Ouillen et alj|2007t) . 

3.1. parameter exploration 

Using {iw) as an approximate measure of the disc's vertical 
thickness (through the relation {iw) ~ ^fi2)h), we display in 
Fig|5]how our results vary when changing the values of the main 
parameters considered here. 



4. Discussion 

4. 1 . limitations and strength of our approach 

We want to stress that this work is in no way an attempt at real- 
istically modelling the coupled dynamical and collisional evo- 
lution of debris discs. As already mentioned, no determinis- 
tic code can presently achieve this task, mainly because a cor- 
rect treatment of collisions would require to take into account 
the numerous individual small fragments produced after each 
high-velocity impact, thus leading to an exponential increase 
of simulated test particles, rapidly exceeding any reasonable 
value 0. Our study should be regarded as a numerical experi- 
ment aimed at exploring one specific mechanism: the ability of 
mutual collisions to convert the high in-plane velocities due to 
small, radiation-pressure aff'ected particles into vertical ones. 

The main and unavoidable problem with our approach is 
that it considers that particles acquire their dynamics through 
a succession of inelastic rebounds. This is obviously unrealis- 
tic, especially for the smallest grains close to the blow-out size 
Scut- In our runs, such grains evolve mostly through rebounds 
with grains their own size (which dominate the total population), 
whereas in reality they should be fragments produced from col- 
lisions on bigger parent bodies (while collisions between ~ Scut 
grains should produce small << Scur debris quickly removed by 
radiation pressure). However, for the specific issue addressed 
here, what matters is only if this departure from reality intro- 
duces a systematic error when estimating the final (;). Exactly 
quantifying such a potential bias is clearly beyond the scope of 
this study, but it is possible to get a first idea. Statistical colli- 
sional evolution model (e.g. The bault & Auger eau 2007) show 
that, for extended debris discs with low dynamical excitation 
({e)dy„ < 0.05), grains close to Scm are mostly produced through 



' Note however the pioneering work of iBeauge & Aarseth 
whose code produced 4 fragments after each shattering collision, 
but had to be rest ricted to only 200 initial bodies, or that of 
Leinhard t & Richardson . (2005), allowing the breakup of large rubble 
piles of gravitationally bo und hard spheres, as well as the more recent 
work by Grigorieva et aP l|2007) who combined the statistical and N- 
body approaches by using numerical "super-particles", but were limited 
to short timescales 
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cratering and shattering collisions between ~ Scut projectiles 
and targets in the ~ 5 x i^,,, to ~ 100 x Scui range. In other words, 
the most common way to produce an Scm fragment is by an im- 
pact involving a previous (and ultimately shat tered) ~ Scut pro- 
j ectile . If we now use the simplified estimate of Stern & Colwelfl 
(1 19971) for the velocity of fragments produced after an impact at 
speed AvQ 



, fkeE, 



kin 



M, 



target 



0.5 



(3) 



where Ekin - 0-5(M,arget + Mp,-oj)^v and fke is the fraction of 
kinetic energy not dissipated at impact, we see that vy> ~ ./^'p^Av 
when Mfarget » Mproj. This means that the fragment velocity is 
comparable to the outcome velocity for an inelastic rebound be- 
tween two equal-sized bodies dissipating (1 - fi^,) of the kinetic 
energy at impact. Thus, by tuning in the value of fke, the dynam- 
ical outcome of a collision between two Sdu particles bouncing 
at speed Av can be made to roughly mimic that of a Scut fragment 
produced by a Av impact of an Scut projectile on a larger target. 
This is why the rebound coefficients and et have been given 
values corresponding to the typical energy loss expected from 
high velocity fragmenting impacts (see Sec l2.2b . We are never- 
theless aware that it is impossible to make this N-body code fully 
bias free in this respect. It is in particular likely that it slightly 
overestimates the global Ekin transfer for small ~ Scut grains, 
since a small fraction of them should in reality be produced by 
fragmenting collisions involving projectiles > for which the 
impact kinetic energy is lower than that for ~ projectiles. 

A second issue is that, due to the succession of dissipative 
collisions, the numerical system's total kinetic energy is steadily 
decreasing. In particular, the high initial (e) of the smaller parti- 
cles are progressively damped as they suffer more and more im- 
pacts. This is an artificial behaviour, since in reality the smaller 
grains should not survive several high-Av impacts and should be 
destroyed before having had the time to be dynamically damped, 
while newly produced grains of the same size should progres- 
sively enter the system with high initial eccentricities ("automat- 
ically" acquired when they are collisionally released from their 
parent body). To minimize this numerical bias, simulations are 
stopped before the average number of collisions per body ex- 
ceeds 3. As can be seen from Fig|2l this is enough for the system 
to reach a steady state in terms of (/) for each grain population. 
Nevertheless, an artificial decrease of the system's excitation is 
observed, which leads to underestimating the final (/). 

We thus see that it is impossible to make our N-body ap- 
proach fully bias free. Nevertheless, for the reasons just dis- 
cussed, the two main sources of systematic errors should remain 
relatively limited. Moreover, they should to some extent com- 
pensate each other As a consequence, we are relatively confident 
that, for systems which are sufficiently evolved to have reached 
a steady state where the size distribution and dynamical struc- 
ture no longer evolves due to collisions, the present code gives a 
satisfying first-order estimate. As mentioned in Sect|2] its main 
strength is that it takes into account the main parameters that 



* this is mainly because grains close to the blow out size are the im- 
pactors carrying the highest kinetic ener gy and have thus an incr eased 
excavating/shattering capacity (see The bault & Augereau II2007I . for a 
detailed discussion on these issues) 

' This is of course a very simplified expression, which neglects 
the fact that v;> should in principle depend on fragment size and on 
the way fragments are produced: cratering or shattering of the tar- 
get. Nevertheless, we can ignore these refinements for the order-of- 
magnitude purpose of our discussion 



drive the in-to-out-of-plane kinetic energy transfer: collision ve- 
locity, impact geometry, mass ratio between the impactors and 
energy loss at impact. Note that such a " bouncing bal l s" ap - 
proach has also been recently proposed by iChiang et al.l (l2009l) 
as the best way to address a not-so-different issue, i.e., the col- 
lisional relaxation of grain orbits in the FomaUiaut ring. 

4.2. the "natural" thickness of debris discs 

Despite the limitations inherent to our numerical approach (see 
previous section), we believe our main result to be relatively ro- 
bust, i.e., there is a "natural" thickening of a debris disc due to 
the combined action of radiation pressure and mutual collisions. 
Expressing the disc's thickness in terms of its average aspect ra- 
tio h - H/R, we find that (from Fig|5] and using the relation 
(iw) ~ ^^(2)h ) 



h,„i„ = 0.04 ± 0.02 



(4) 



where the error bar is probably overestimated, as it is obtained 
when assuming rather extreme values for our parameter explo- 
ration (see Sec B.lb . A direct consequence is that there should be 
in principle no debris disc with an aspect ratio smaller than /!,„,„, 
at least in the generic case of systems whose geometrical cross 
section is dominated by the smallest grains close to the blow-out 
size (see discussion ot the end of this subsection). 

Another important result is that the disc's vertical thickness 
is almost independent of its intrinsic dynamical excitation as 
long as (edyn) - 2 < 0.15. In this case, all intrinsic dy- 
namical effects, in particular the out of plane excursion due to 
the "dynamical" /, are masked by the "natural" thickening ef- 
fect studied here. Likewise, it means that for a disc with an ob- 
served h < 0.04 (or 0.06 to be conservative), the coiTelation be- 
tween vertical structure and intrinsic dynamical excitation is lost. 
In particular, no information on potential embedded perturbing 
objects can be retrieved, except that they cannot be larger than 
~ 500 km (EqlB. 

A third result is that discs should be stratified with respect 
to sizes, with the bigger grains being clustered much closer to 
the midplane. This might have interesting observational conse- 
quences. Indeed, it should imply that discs should appear much 
thinner at longer wavelengths, where smaller grains are poor 
scatterer/emitter and where the flux is dominated by bigger par- 
ticles. For a typical value Scut ~ 2 - 10/im, this should thus hap- 
pen for A > 50/im, i.e. in the mid to far IR thermal emission. 
Unfortunately high-resolution images are very difficult to obtain 
at these long wavelengths and there is to our knowledge no de- 
bris disc that is vertically resolved beyond the mid-IR. 

There might however be one marginal possibility for a disc to 
appear flatter than /i„„„ even in scattered light, i.e. if its total geo- 
metrical cross section is dominated by grains significantly larger 
than the blow out size. This is clearly not the case for a Dohnanyi 
size distribution in dN cc s^^'^ds, as well as for most of the more 
realistic "wa vy" distributions derived from coUisional evolutio n 
models (e.g. iKrivov et al.1 l2006t iThebault & Augereau I l2007l) . 
There is however a range of systems for which these numerical 
studies have shown that they could be depleted in small grains, 
i.e., those wi th very low dynamical excitation (edyn)- Indeed, as 
identified bv IThebault & Wul (l2008l) . for (edy„) < 0.01, there is 
an imbalance between the rates of small grain production, which 
strongly decreases with edvri, and destruction, which only mod- 
erately changes for low e^y„ values (because it is mostly imposed 
by spr). For these systems, the scattered light flux should be 
dominated by grains typically 10-100ie,„ in size, which should 
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Fig. 6. Vertical F WHM (^2H) for th e AU Mic disc, as derived 
from observations (iKrist et al.ll2005l) and obtained for a test nu- 
merical disc initially confined to a narrow annulus at ~ 40 AU 
with (erfvn) = (idyn) = 0. The observed FWHM has bee n av- 
eraged over the two ansae (see Fig.7 of iKrist et aUlIOOl) . The 
synthetic profile has been degraded to the resolution of the scat- 
tered light images, i.e. 2 AU 



have very lo w ; and thus lead to a very flat aspect. Note how- 
ever that the iThebault & Wul (|2008|) estimates did not take into 
account the mechanism identified here of collisional redistribu- 
tion of in-plane to out-of -plane velocities, which should tend to 
reduce the rate of small grains destructions by spatially diluting 
them in z- It remains to see to which extent this effect might re- 
duce the depletion of ~ i^. ,„ grains and if these dynamically cold, 
small-grains-depleted systems are a realistic alternative. 

Let us stress that this discussion regards the global aspect 
ratio or thickness of the disc, not possible local clumps, blobs 
or other asymmetries in its vertical structure, like the famous 
beta Pic warp. Such local features could well be the signature of 
hidden perturbers but are not the purpose of the present work. 

4.3. comparison to real systems 

As already mentioned, the two best resolved debris discs, 
y6 Pictoris and AU Mic0, have aspect ratios h ~ 0.05-0. 1 and 0.05 
respectively. AU Mic's vertical thickness thus falls within the 
range for the "natural" /i„„„ derived in EqH) whereas /3 Pictoris 
is around the upper end of this range. Note that the preliminary 
estimate for HD 15 115 aspect's ratio, i.e. ~ 0.05, also falls within 
our estimated range for /?,„,„. 

To push our analysis a step further we perform one additional 
run specifically aimed at fitting AU Mic (/? Pictoris has a much 
more complex radial and vertical structure, especially its warp, 
which does not make it an ideal candidate for such a simple fit). 
For this task, we take input parameters adapted to the specifici- 
ties of this system, i.e., a "birth ring" (where all non radiation 
affected larger par ticles are loca ted) centered at ~ 40 AU and 
of width ~ 5 AU (lAugereau & B eust 2006; S trubbe & Chiand 



l2006l) . As for our previous runs, we consider a limiting dynami- 
cally cold case with (erfv„) = (/,■„,>,«/) = 0, the rest of the param- 
eters are those for our nominal set-up (TablT]i. As can be seen 
in Fig|6] a reasonably good fit of the system's FWHM radial 
profile is obtained. Note that this simple fit has been obtained 
with our arbitrarily defined nominal set-up, and we expect much 
better fits to be possible when tuning in all free parameters (re- 
bound coefficients, size distribution, etc.). Of course, this result 
should not be regarded as a proof that the AU Mic disc is com- 
pletely dynamically cold. It is just a numerical test showing that 
this system's vertical structure is easily explained by the natural 
thickening mechanism identified in this study, without the need 
of additional exciting sources such as hidden large planetesimals 
or planets. The presence of such hidden objects can of course not 
be ruled out, and they might even make sense within the frame of 
our current understanding of debris discs, but we want to stress 
that nothing can be inferred about their presence or absence from 
the disc's observed aspect ratio. 

There is however one specific system which might appear as 
a counter example to our study: FomaUiaut and its striking debris 
ring surrounding a recently discovered giant planet (Kalas et aT| 
2008). Indeed, this ring's aspect ratio has been observ ationally 
estimated to be /z = H/r ~ 0.0125 ( Kalas et al.]|2005l) . a value 
which was corrobo rated by the numerical estimate h ~ 0.017 by 
IChiang et all (12009) . These values are below the minimum value 
for hmin derived in Equ|4]and seem to invalidate our conclusions. 
However, two important points should be stressed here. The first 
one is that the disc's thickness estimates were obtained for one 
specific region, i.e., the sharp inner edge of the bright ring lo- 
cated at 133 AU0. The second one is that the thickness which 
has been estimated for Fomalhaut is the local geometrical verti- 
cal thickness Hg^om, and not the line of sight integrated thickness 
of the luminosity profile as for an edge-on seen disc. We have 
seen in Sect|3]that, although these two different thicknesses do 
most of the time give similar values, Hge„,n distinguishes itself 
in one specific case, i.e., the inner edge of the parent bodies disc 
where it can reach very small values (see FiglDl. Exactly how 
small is difficult to evaluate because of the limited radial reso- 
lution of our synthetic profiles, but it is clearly below //,„,„, and 
could thus be compatible with the Fomalhaut estimate. Clearly, 
this specific case has to be investigated in more details. It would 
in particular be of great help to obtain a thickness estimate fur- 
ther out in the disc, where our model predicts a wider vertical 
dispersion. 



5. Summary and Conclusions 

The main results of our numerical investigation can be summa- 
rized as follows: 

- Even in the absence of perturbing bodies, debris disc should 
naturally thicken because of the coupled effect of radiation 
pressure and mutual grain collisions 

- The basic mechanism is the following: because of radiation 
pressure, the smallest grains have high in-plane velocities 
which are partially converted into vertical motions by col- 
lisions with other grains. The dynamical outcome depends 
on the grain differential sizes and the collision geometry. 



AU Mic is an M star where, in addition to radiation pressure, stellar 
wind is also acting to place small grains on e c centric or unbound orbits 
dWood et al.ll200i lAugereau & Beuslll2006l : IStrubbe & Chiang! l2006l : 
[Fitzgerald et al. 200"^. This additional force has to a first order the same 
dependency on particle size and radial distance as radiation pressure. 



' They are based on the measure of this edge's drop-off in surface 
brightness, which, because of the system's inclination to the line of 
sight, could either be due to an infinitely flat disc with a finite spread in 
the radial direction, or a per fectly sharp ra dial edge with a finite vertical 
thickness (see Sect.3.2.2 of IChiang et"aL|[2"009l) 
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- The disc is stratified in the vertical direction with respect to 
particle sizes, the smallest grains having the largest disper- 
sion while the bigger ones remain close to the midplane. 

- At wavelengths were the flux is dominated by the smallest 
grains, i.e., in the visible to mid-IR, this natural thickening 
effect places a minimum value /i,„,„ on a disc's aspect ratio. 
We find h„,i„ ~ 0.04 ± 0.02, although this estimate should 
be taken with some care given the unavoidable limitations of 
our model. 

- Most vertically resolved debris discs have aspect ratios com- 
patible with h,„i„, which can thus be explained without re- 
sorting to the perturbing presence of large "hidden" bodies. 
As an illustration, we obtain a relatively satisfying fit of AU 
Mic's vertical profile. 

- A corollary is that, for all systems with h ~ h,„in, no infor- 
mation about the presence of embedded planetary embryos 
can be directly retrieved from measuring the disc's vertical 
dispersion. 

Let us conclude by stressing that these results should in no way 
be interpreted as a proof that debris discs are dynamically cold 
and devoid of hidden massive perturbing bodies. Indeed, the 
presence of such perturbers is fully compatible with our results. 
As an example, a disc with h - 0.04 could harbour embryos as 
large as a few 100 km. As a matter of fact, we subscribe to the 
theoretical arguments supporting the presence of large hidden 
bodies, which are in particular needed to provide the mass reser- 
yoir for sustaining the high collisional activity of debris discs 
jLohne et al.l 120081) . What we have shown is that direct obser- 
vational evidence can probably not be easily obtained, since for 
aspect ratios lower than ~ 0.06 there is no direct link between 
vertical thickness and the dynamical excitation of the system. 
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List of Objects 

'ySPictoris' on page 7 
'jSPictoris' on page 7 
'jSPictoris' on pageQ 



